Prescribed burns maps

# Loading required packages for analysis:
library(here)
library(tidyverse)
library(ggplot2)
library(sf)
library(ggspatial)
library(gganimate)
library(gifski)
library(transformr)
library(ggrepel)

source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
source(here("Scripts", "4_Covariate_Data_Prep.R"))
# source(here("Scripts", "5_Prep_Species.R"))

# Generate basemap (habitat types)
hab_groups_forest <- c("Oak Woods", "Early Successional Forest", "Shrubland", "Maple Basswood", "Flood Plain Forest", "Tamarack Bog")
hab_groups_grasslands <- c("Oldfield", "Prairie")
hab_groups_wetlands <- c("Wetland Cattail", "Wetland Sedge", "Wetland Lily", "Beach Mudflat")
hab_group_water <- c("Open Water")
hab_groups_cultivated <- c("Planted Cultivated")

basemap_habs <- spatial.habitats.current %>%
  mutate(type_simplified = ifelse(MLCCS_Data %in% hab_groups_forest, "forest",
    ifelse(MLCCS_Data %in% hab_groups_grasslands, "grassland",
      ifelse(MLCCS_Data %in% hab_groups_wetlands, "wetland",
        ifelse(MLCCS_Data %in% hab_groups_cultivated, "cultivated",
          ifelse(MLCCS_Data %in% hab_group_water, "water", NA)
        )
      )
    )
  )) %>%
  filter(!is.na(type_simplified))

basemap_cols <- c("forest" = "darkseagreen", "grassland" = "lightgoldenrod", "cultivated" = "wheat3", "wetland" = "lightblue", "water" = "skyblue4")

basemap <- ggplot() +
  geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE) +
  scale_fill_manual(values = basemap_cols) +
  theme_void()

basemap

# Generate park name lookup table for parks with burns (relates codes in burn data to park names in park park boundaries dataset).
park_name_lookup <- data.frame(
  park = c("BAK", "BRY", "CAR", "CRO", "ELM", "FR", "GW", "HYL", "KW", "MUH", "MUL", "NHTC", "REB", "SPR"),
  Name = c("Baker", "Bryant", "Carver", "Crow-Hassan", "Elm Creek", "French", "Gale Woods", "Hyland", "Kingswood", "Murphy-Hanrehan", "The Landing", "NHTC", "Lake Rebecca", "Spring Lake Park")
)

# Function for quickly defining map extent in ggplot based on a particular park. Park name must match that listed in park boundaries sf.
map_extent_park <- function(park) {
  row <- which(spatial.boundaries.parks$Name == park)
  bbox <- st_bbox(spatial.boundaries.parks[row, ])
  coord_sf(
    xlim = c(bbox$xmin, bbox$xmax),
    ylim = c(bbox$ymin, bbox$ymax)
  )
}

# Function for mapping burns in a particular park/year
map_burns <- function(park, display_year) {
  units_filter_year <- filter(units, year == display_year)
  cols <- c("1" = "orangered3", "forest" = "darkseagreen", "grassland" = "lightgoldenrod", "wetland" = "lightblue", "water" = "skyblue4", cultivated = "wheat3")
  ggplot() +
    geom_sf(data = spatial.boundaries.parks, fill = "grey97", size = 0) +
    geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE, alpha = 0.7) +
    # geom_sf(data = spatial.infrastructure.trails %>%
    #               filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
    #               st_zm(drop = TRUE),
    #        linetype = "dotted", size = .5, alpha = .8, color = "tan4") +
    geom_sf(
      data = spatial.infrastructure.roads %>%
        filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
        st_zm(drop = TRUE),
      linetype = "solid", size = .4, alpha = .8, color = "tan4"
    ) +
    geom_sf(data = units_filter_year %>% filter(burn_status_simple == 1), aes(fill = as.factor(burn_status_simple)), show.legend = FALSE, color = "white", alpha = 0.75) +
    geom_sf(data = spatial.boundaries.parks, fill = NA, size = 0.7) +
    scale_fill_manual(values = cols) +
    scale_y_continuous(breaks = seq(40, 50, by = 0.02)) +
    scale_x_continuous(breaks = seq(-95, -90, by = 0.02)) +
    theme_minimal() +
    theme(
      axis.title = element_blank(),
      plot.title = element_text(size = 11, face = "bold", color = "orangered3"),
      plot.subtitle = element_text(size = 10, face = "plain")
    ) +
    map_extent_park(park) +
    annotation_scale(location = "br", style = "bar", bar_cols = c("grey40", "grey40"), line_col = "grey40", text_col = "grey40", line_width = 0.5) +
    labs(
      title = "Prescribed burns",
      subtitle = paste0(park, ": ", display_year)
    )
}

map_burns("Crow-Hassan", "1990")

maps_2023 <- lapply(park_name_lookup$Name, map_burns, "2023")
names(maps_2023) <- park_name_lookup$Name
maps_2023

# Function for generating an animation of burns in a particular park
map_burns_for_ani <- function(park) {
  # Define colors
  cols <- c("1" = "orangered3", "cultivated" = "wheat3", "forest" = "darkseagreen", "grassland" = "lightgoldenrod", "wetland" = "lightblue", "water" = "skyblue4")

  # Pull burned untis only
  units_burned <- units %>% filter(burn_status_simple == 1)

  # Get park code from park name
  park_code <- park_name_lookup$park[which(park_name_lookup$Name == park)]

  # Calculate the minimum start year for the current park
  min_start_year <- min(as.numeric(as.character(units_burned$year[which(units_burned$park == park_code)])))

  # Create a subset of data starting from the year before the first burn
  units_subset <- units_burned %>% filter(park == park_code, as.numeric(as.character(year)) >= min_start_year - 1)

  # For mapping purposes, we need a fake/placeholder row in the birn units database (one for each unique combination of park and year). I'll have the geometry be a small circle centered on Minneapolis (out of frame of any park).
  fake_geometry <- st_sfc(st_point(c(-93.2650, 44.9778)))
  fake_geometry <- st_buffer(fake_geometry, dist = 100)
  st_crs(fake_geometry) <- st_crs(units)

  # Create a data frame with unique combinations of year and park
  fake_data <- expand.grid(year = unique(units$year), park = unique(units$park), burn_status_simple = 1, unit = NA, Burn_unit_description = NA) %>%
    filter(as.numeric(as.character(year)) >= min_start_year - 1) %>%
    filter(park == park_code)

  # Combine the data frame with the spatial data
  fake_df <- st_sf(fake_data, Shape = fake_geometry)

  # Bind the fake rows to the original units_burned dataframe
  units_subset <- rbind(units_subset, fake_df)

  # Generate basemap for animation
  plot <- ggplot() +
    geom_sf(data = spatial.boundaries.parks, fill = "grey97", size = 0) +
    geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE, alpha = 0.7) +
    # geom_sf(data = spatial.infrastructure.trails %>%
    #               filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
    #               st_zm(drop = TRUE),
    #        linetype = "dotted", size = .5, alpha = .8, color = "tan4") +
    geom_sf(
      data = spatial.infrastructure.roads %>%
        filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
        st_zm(drop = TRUE),
      linetype = "solid", size = .4, alpha = .8, color = "tan4"
    ) +
    geom_sf(data = units_subset, aes(fill = as.factor(burn_status_simple)), show.legend = FALSE, color = "white", alpha = 0.75) +
    geom_sf(data = spatial.boundaries.parks, fill = NA, size = 0.7) +
    scale_fill_manual(values = cols) +
    scale_y_continuous(breaks = seq(40, 50, by = 0.02)) +
    scale_x_continuous(breaks = seq(-95, -90, by = 0.02)) +
    theme_minimal() +
    theme(
      axis.title = element_blank(),
      plot.title = element_text(size = 11, face = "bold", color = "orangered3"),
      plot.subtitle = element_text(size = 10, face = "plain")
    ) +
    map_extent_park(park) +
    annotation_scale(location = "br", style = "bar", bar_cols = c("grey40", "grey40"), line_col = "grey40", text_col = "grey40", line_width = 0.5) +
    labs(
      title = "Prescribed burns",
      subtitle = paste0(park, ": {next_state}")
    ) # not sure why next_state works but closest_State doesn't; closest lags the data by 1 year

  # Facet by year for animation frames
  plot_animation <- plot + transition_states(
    year,
    transition_length = 0,
    state_length = 1,
    wrap = TRUE
  )

  # Generate animation
  mapGIF <- animate(
    plot_animation,
    nframes = length(unique(units_subset$year)) + 2,
    fps = 1.5,
    end_pause = 2,
    renderer = gifski_renderer()
  )

  # Save animation
  anim_save(
    filename = paste0("burns_", park, ".gif"),
    animation = mapGIF,
    path = here("Results", "990_burn_animation")
  )
}

parks_list <- c("Baker", "Bryant", "Carver", "Crow-Hassan", "Elm Creek", "French", "Gale Woods", "Hyland", "Kingswood", "Murphy-Hanrehan", "The Landing", "Lake Rebecca", "Spring Lake Park") # parks with  burn data
lapply(parks_list, map_burns_for_ani)

files <- list.files(path = here("Results", "990_burn_animation"), pattern = ".gif", full.names = TRUE)
files_names <- list.files(path = here("Results", "990_burn_animation"), pattern = ".gif", full.names = FALSE)
files_names <- gsub("\\.gif", "", files_names)
files_names <- gsub("burns_", "", files_names)

for (i in 1:length(files)) {
  cat(paste0("#### ", files_names[i], "\n"))
  cat(paste0("![](", files[i], ")\n\n"))
}